library(tidyverse)
library(dplyr)
library(tibble)
library(haven)
library(sandwich)
library(MASS)
library(ebal)
library(MatchIt)
library(cobalt)
library(lmtest)
library(Matching)
library(effects)
library(cat)
library(splines)
library(jtools)
library(broom)
library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(pastecs)
library(stargazer)
library(WeightIt)
library(modelsummary)
library(margins)

## Read data 
merge1 <- read.csv("MainData.csv")
raw_capex <- read.csv("raw_capex.csv")
winsorized <- read.csv("winsorized.csv")

## Main models 
set.seed(1000)
#Full matching 
fm1 <- matchit(treatment ~ US_allies + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia,
               data = merge1, method = "full", 
               distance = "glm", link = "probit")
summary(fm1)
fmm1 <- match.data(fm1)
factor(fmm1$Grade, order = TRUE, 
       levels = c("Digging In", "Buying Time", "Scaling Back", "Suspension", "Withdrawal")) -> fmm1$Grade


#Factorize for Industry FEs and for plotting
fmm1$Industry_Sector <- as.factor(fmm1$Industry_Sector)
fmm1$extractive_firm <- as.factor(fmm1$extractive_firm)
bal.tab(fm1)

#Proportions of firms in each category (Figure 1)
ggplot(data = merge1, aes(x = as.factor(Grade), fill = Industry_Sector)) +
  geom_bar(position = "fill") +
  labs(x = "Grade", y = "Proportion") +
  scale_y_continuous(labels = scales::percent)

#Common Support (Figure 2)
fmm1 %>% 
  mutate(Treatment = ifelse(treatment == 1, "Treatment", "Control")) %>% 
  ggplot()+
  geom_histogram(aes(x = distance, y=..density.., fill= Treatment), color = "black", bins = 30, position = "identity", alpha=.3) +
  geom_density(aes(x = distance, color = Treatment), alpha = .6) +
  ggtitle("Common Support") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) 


#Baseline models (Table 1)
m1_bin <- glm(Withdraw ~ Russia + US_allies + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia + Industry_Sector, data = fmm1, weights = weights, family = "binomial")
m1_bin_cl <- coeftest(m1_bin, cluster = ~  Country + subclass)
print(m1_bin_cl)

m2_bin <- glm(Withdraw ~ extractive_firm + US_allies + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia, data = fmm1, weights = weights, family = "binomial")
m2_bin_cl <- coeftest(m2_bin, cluster = ~  Country + subclass)
print(m2_bin_cl)
stargazer(m1_bin_cl, m2_bin_cl, omit = c("Industry_Sector"))

#Economic footprint in relative term (with winsorization to account for outliers)
fmm1$relative_footprint <- winsorized$x
m1_relative <- glm(Withdraw ~ relative_footprint + US_allies + ln_employee + ln_market_cap + year_in_Russia + Industry_Sector, data = fmm1, weights = weights, family = "binomial")
m1_relative_cl <- coeftest(m1_relative, cluster = ~  Country + subclass)
print(m1_relative_cl)

#Plotting 
#Figure 3 (Balancing Plot)
love.plot(fm1, abs = TRUE, drop.distance = TRUE, var.order = "unadjusted", stars = "raw",
          line = FALSE)

#Figure 4 (H1)
plot_model(m1_bin, type = "pred", terms = "Russia [0:66]") + ggtitle("Marginal Effects of Russia Subsidiaries on Withdrawal Likelihood") + 
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) + 
  xlab("Russia_Subsidiaries")

#Figure 5 (H2)
plot_model(m2_bin, type = "pred", terms = "extractive_firm") + ggtitle("Predicted Probability of Withdrawal") + 
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  xlab("extractive")


## Robustness test/Appendix 
#List of extractive-related firms (Table 1)
merge1 %>% 
  filter(extractive_firm == 1) %>% 
  dplyr::select(Name, Industry_Sector) -> list_of_firms

latex_code <- xtable::xtable(list_of_firms)
print(latex_code)

#Summary Statistic (Table 2)
merge1 %>% 
  dplyr::select(extractive_firm, Russia, Sub_US_allies, US_allies, ln_market_cap, ln_employee, year_in_Russia) -> sum_table
stargazer(sum_table)

#Balancing Equation (Table 4)
bal_eq <- glm(treatment ~ US_allies + ln_employee + Sub_US_allies + year_in_Russia + ln_market_cap, data = merge1, family = "binomial")
summary(bal_eq)
stargazer(bal_eq)

#Entropy Balancing (Table 5)
bal_model1 <- ebalance(Treatment = merge1$treatment, X = merge1[, c("US_allies", "Sub_US_allies", "ln_employee", "year_in_Russia", "ln_market_cap")])
control1 <- merge1 %>% 
  filter(treatment == 0) %>% 
  mutate(weights = bal_model1$w)
treated1 <- merge1 %>% 
  filter(treatment == 1) %>% 
  mutate(weights = 1)
df_bal1 <- bind_rows(control1, treated1)

df_bal1$extractive_firm <- as.factor(df_bal1$extractive_firm)
m1_bal <-  glm(Withdraw ~ Russia + US_allies + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia + Industry_Sector, data = df_bal1, weights = weights, family = "binomial")
m1_bal_cl <- coeftest(m1_bal, cluster = ~Country)
m2_bal <- glm(Withdraw ~ extractive_firm + US_allies + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia, data = df_bal1, weights = weights, family = "binomial")
m2_bal_cl <- coeftest(m2_bal, cluster = ~ Country)

#BIT (Table 5)
m1_bin_BIT_cl <- coeftest(glm(Withdraw ~ Russia + US_allies + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia + Industry_Sector + BIT, data = fmm1, weights = weights, family = "binomial"), cluster = ~  Country + subclass)
m2_bin_BIT_cl <- coeftest(glm(Withdraw ~ extractive_firm + US_allies + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia + BIT, data = fmm1, weights = weights, family = "binomial"), cluster = ~  Country + subclass)

stargazer(m1_bal_cl, m2_bal_cl, m1_bin_BIT_cl, m2_bin_BIT_cl, omit = c("Industry_Sector"))

#Ordered logit (Table 6)
factor(fmm1$Grade, order = TRUE, 
       levels = c("Digging In", "Buying Time", "Scaling Back", "Suspension", "Withdrawal")) -> fmm1$Grade

m1_1 <- polr(Grade ~ Russia + US_allies + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia + Industry_Sector, data = fmm1, weights = weights, method = "logistic")
m1_1_cl <- coeftest(m1_1, cluster = ~ Country + subclass) 

ol <-  polr(Grade ~ extractive_firm + US_allies + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia, data = fmm1, weights = weights, method = "logistic")
ol_cl <- coeftest(ol, cluster = ~ Country + subclass)
stargazer(m1_1_cl, ol_cl, omit = c("Industry_Sector"))

#Full matching on US_firm (Table 7)
fm2 <- matchit(treatment ~ US_firm + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia,
               data = merge1, method = "full", 
               distance = "glm", link = "probit")
fmm2 <- match.data(fm2)

m1_fmm2 <- coeftest(glm(Withdraw ~ Russia + US_firm + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia + Industry_Sector, data = fmm2, weights = weights, family = "binomial"), cluster = ~ Country + subclass)
m2_fmm2 <- coeftest(glm(Withdraw ~ extractive_firm + US_firm + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia, data = fmm2, weights = weights, family = "binomial"), cluster = ~ Country + subclass)
stargazer(m1_fmm2, m2_fmm2, omit = c("Industry_Sector"))

#SOEs (Table 8)
SOE <- read.csv("SOE.csv")

m3_bin <- glm(Withdraw ~ Russia + US_allies + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia + Industry_Sector + SOE, data = SOE, weights = weights, family = "binomial")
m3_bin_cl <- coeftest(m3_bin, cluster = ~  Country + subclass)

m3_bin_extractive <- glm(Withdraw ~ extractive_firm + US_allies + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia + SOE, data = SOE, weights = weights, family = "binomial")
m3_bin_extractive_cl <- coeftest(m3_bin_extractive, cluster = ~  Country + subclass)

stargazer(m3_bin_cl, m3_bin_extractive_cl, omit = c("Industry_Sector"))


#Baseline models without matching (Table 9)
m_simple <- glm(Withdraw ~ Russia + US_allies + ln_employee + Sub_US_allies + ln_market_cap + Industry_Sector + year_in_Russia, data = merge1, family = "binomial")
m_simple_cl <- coeftest(m_simple, cluster = ~ Country)
print(m_simple_cl)

m_simple1 <- glm(Withdraw ~ extractive_firm + US_allies + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia, data = merge1, family = "binomial")
m_simple1_cl <- coeftest(m_simple1, cluster = ~ Country)
print(m_simple1_cl)
stargazer(m_simple_cl, m_simple1_cl, omit = c("Industry_Sector"))

#Ordered logit without matching (Table 10)
factor(merge1$Grade, order = TRUE, 
       levels = c("Digging In", "Buying Time", "Scaling Back", "Suspension", "Withdrawal")) -> merge1$Grade 
m_simple_ordered <- polr(Grade ~ Russia + US_allies + ln_employee + Sub_US_allies + ln_market_cap + Industry_Sector + year_in_Russia, data = merge1, method = "logistic")
m_simple_ordered_cl <- coeftest(m_simple_ordered, cluster = ~ Country)
m_simple_ordered1 <- polr(Grade ~ extractive_firm + US_allies + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia, data = merge1, method = "logistic")
m_simple_ordered1_cl <- coeftest(m_simple_ordered1, cluster = ~ Country)

stargazer(m_simple_ordered_cl, m_simple_ordered1_cl, omit = c("Industry_Sector"))

#Alternative measurement for economic footprint (Table 11)
fmm1 %>% 
  filter(Name %in% raw_capex$Name) -> subset_fmm1 

m_capex <- glm(Withdraw ~ ln_capex + US_allies + ln_employee + Sub_US_allies + ln_market_cap + year_in_Russia + factor(Industry_Sector), data = subset_fmm1, weights = weights, family = "binomial")
m_capex_cl <- coeftest(m_capex, cluster = ~ Country + subclass)
stargazer(m_capex_cl, omit = c("Industry_Sector"))

#Interaction model (Table 12)
fmm1$after_2014 <- ifelse(fmm1$Year > 2014, 1, 0)
m_int <- glm(Withdraw ~ Russia*after_2014 + US_allies + ln_employee + Sub_US_allies + ln_market_cap + Industry_Sector, data = fmm1, weights = weights, family = "binomial")
m_int_cl <- coeftest(m_int, cluster = ~  Country + subclass)
stargazer(m_int_cl, omit = c("Industry_Sector"))

